NATIONAL  BUREAU  OF  STANDARDS  REPORT 


A STATES TECAL  METHOD  FOR  DETERMINING  THE  LOWEST 


EIGENVALUE  OF  SCHRODINGER 1 S EQUATION 


by 

Mark  Kac 

Cornell  University 
and 

National  Bureau  of  Standards 
and 

Michael  Cohen 

Cornell  University 
and 

National  Bureau  of  Standards 


U.  S.  DEPARTMENT  OF  COMMERCE 
NATIONAL  BUREAU  OF  STANDARDS 


1#3 


U.  S.  DEPARTMENT  OF  COMMERCE 
Charles  Sawyer,  Secretary 


NATIONAL  BUREAU  OF  STANDARDS 
A.  V.  Astin,  Director 

THE  NATIONAL  BUREAU  OF  STANDARDS 

The  scope  of  activities  of  the  National  Bureau  of  Standards  is  suggested  in  the  following  listing  of 
the  divisions  and  sections  engaged  in  technical  work.  In  general,  each  section  is  engaged  in  special- 
ized research,  development,  and  engineering  in  the  field  indicated  by  its  title.  A brief  description 
of  the  activities,  and  of  the  resultant  reports  and  publications,  appears  on  the  inside  of  the  back 
cover  of  this  report. 

Electricity.  Resistance  Measurements.  Inductance  and  Capacitance.  Electrical  Instruments. 
Magnetic  Measurements.  Applied  Electricity.  Electrochemistry. 

Optics  and  Metrology.  Photometry  and  Colorimetry.  Optical  Instruments.  Photographic 
Technology.  Length.  Gage. 

Heat  and  Power.  Temperature  Measurements.  Thermodynamics.  Cryogenics.  Engines  and 
Lubrication.  Engine  Fuels.  Cryogenic  Engineering. 

Atomic  and  Radiation  Physics.  Spectroscopy.  Radiometry.  Mass  Spectrometry.  Solid  State 
Physics.  Electron  Physics.  Atomic  Physics.  Neutron  Measurements.  Infrared  Spectroscopy. 
Nuclear  Physics.  Radioactivity.  X-Rays.  Betatron.  Nucleonic  Instrumentation.  Radio- 
logical Equipment.  Atomic  Energy  Commission  Instruments  Branch. 

Chemistry.  Organic  Coatings.  Surface  Chemistry.  Organic  Chemistry.  Analytical  Chemistry. 
Inorganic  Chemistry.  Electrodeposition.  Gas  Chemistry.  Physical  Chemistry.  Thermo- 
chemistry. Spectrochemistry.  Pure  Substances. 

Mechanics.  Sound.  Mechanical  Instruments.  Aerodynamics.  Engineering  Mechanics.  Hy- 
draulics. Mass.  Capacity,  Density,  and  Fluid  Meters. 

Organic  and  Fibrous  Materials.  Rubber.  Textiles.  Paper.  Leather.  Testing  and  Specifi- 
cations. Polymer  Structure.  Organic  Plastics.  Dental  Research. 

Metallurgy.  Thermal  Metallurgy.  Chemical  Metallurgy.  Mechanical  Metallurgy.  Corrosion. 

Mineral  Products.  Porcelain  and  Pottery.  Glass.  Refractories.  Enameled  Metals.  Con- 
creting Materials.  Constitution  and  Microstructure.  Chemistry  of  Mineral  Products. 

Building  Technology.  Structural  Engineering.  Fire  Protection.  Heating  and  Air  Condition- 
ing. Floor,  Roof,  and  Wall  Coverings.  Codes  and  Specifications. 

Applied  Mathematics.  Numerical  Analysis.  Computation.  Statistical  Engineering.  Machine 
Development. 

Electronics.  Engineering  Electronics.  Electron  Tubes.  Electronic  Computers.  Electronic 
Instrumentation. 

Radio  Propagation.  Upper  Atmosphere  Research.  Ionospheric  Research.  Regular  Propaga- 
tion Services.  Frequency  Utilization  Research.  Tropospheric  Propagation  Research.  High 
Frequency  Standards.  Microwave  Standards. 

Ordnance  Development.  These  three  divisions  are  engaged  in  a broad  program  of  research 

Electromechanical  Ordnance,  and  development  in  advanced  ordnance.  Activities  include 
Electronic  Ordnance.  basic  and  applied  research,  engineering,  pilot  production,  field 

testing,  and  evaluation  of  a wide  variety  of  ordnance  materiel.  Special  skills  and  facilities  of  other 
NBS  divisions  also  contribute  to  this  program.  The  activity  is  sponsored  by  the  Department  of 
Defense. 

Missile  Development.  Missile  research  and  development:  engineering,  dynamics,  intelligence, 
instrumentation,  evaluation.  Combustion  in  jet  engines.  These  activities  are  sponsored  by  the 
Department  of  Defense. 

• Office  of  Basic  Instrumentation  • Office  of  Weights  and  Measures. 


NATIONAL  BUREAU  OF  STANDARDS  REPORT 


NBS  PROJECT 


NBS  REPORT 


1101-11-5100 


March  17,  1952 


1553 


A STATISTICAL  METHOD  FOR  DETERMINING  THE  LOWEST 
EIGENVALUE  OF  SCHRODINGER ' S EQUATION* 

tflT 

Mark  Kac 

Cornell  University 
• and 

National  Bureau  of  Standards 


and 

Michael  Cohen 
Cornell  University 
and 


National  Bureau  of  Standards 


PREPRINT 


*The  preparation  of  this  paper  was  sponsored  (in  part)  by  the 

Office  of  Naval  Research 


The  publication,  reprlr 
unless  permission  is  ob 
25,  D,  C,  Such  permls 
cally  prepared  If  that 


Approved  for  public  release  by  the 
Director  of  the  National  Institute  of 
Standards  and  Technology  (NIST) 
on  October  9,  2015 


part,  Is  prohibited 
ndards,  Washington 
ort  has  been  speclfi- 
>rt  for  Its  own  use. 


Table  of  Contents 


Page 

I . Summary  1 

II.  Heuristic  outline  of  the  theory  1 

III.  Estimation  of  the  error  involved  in  the  method 9 

IV.  Extension  to  the  hydrogen  atom * • 16 

V.  The  slope  method;  results  for  hydrogen,  helium,  and 

n double  hydrogen”  , . 22 

VI.  Techniques  for  speeding  the  sampling  process;  high  speed 

computing  and  importance  sampling 3k 

Appendices  .......  ••  

As  Computation  of  discretization  error  for  the  harmonic 

oscillator U6 


Bs  Further  statistical  properties  of  the  random  walks  * * * 52 
C«  Statistical  method  for  finding  the  lowest  eigenfunction  . 58 
Bibliography 


61 


. 

- ...  ... 

£ • ; I - 


OQ  ! 


NBS  1553 

3-17-52 

PREPRINT  COPI 


A Statistical  Method  for  Detennining  the  Lowest 

••  ^ 

Eigenvalue  of  Schrodinger!s  Equation 

by 

Mark  Kac 

Cornell  University 
and 

National  Bureau  of  Standards 
and 

Michael  Cohen 

Cornell  University 
and 

National  Bureau  of  Standards 


I.  SUMMARY 

The  experiments  discussed  here  are  a continuation  of  those  reported  by 
Donsker  and  Kac  in  [2],  Atomic  potentials,  mainly  that  of  hydrogen,  are 
dealt  with  here.  Improved  techniques  for  accumulating  and  treating  the  data 
are  presented;  coupled  with  a high  speed  digital  computer,  these  techniques 
should  make  possible  the  rapid  determination,  to  within  a few  percent,  of 
the  lowest  eigenvalue  of  many  potentials#  Considerable  attention  is  given 
to  the  errors  in  the  method,  and  an  elementary  account  of  the  underlying 
theory  is  included. 

II.  HEURISTIC  OUTLINE  OF  THE  THEORT 
With  the  introduction  of  an  imaginary  time  variable,  the  time-depen- 
dent Schrodinger  equation  for  a single  particle  in  a one-dimensional 

The  preparation  of  this  paper  was  sponsored  (in  part)  by  the  Office 
of  Naval  Research. 

Now  at  the  California  Institute  of  Technology 
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potential  field  V(x)  can  be  brought  into  the  form 

(1)  — = 1 2 - v(x)P(x,t)  . 

3 1 ^ Dtl 

We  shall  shortly  see  that  this  form  of  the  equation  admits  of  a simple 
probabilistic  interpretation.  There  is  some  instructive  value  in  the  pro- 
bability approach,  inasmuch  as  it  adds  to  our  feeling  for  the  meaning  of 
the  wave-function  and  the  equation.  Furthermore,  the  probability  model  for 

equation  (l)  yields  a new  computational  method  for  finding  the  lowest  eigen* 

»• 

value  of  the  time-independent  Schrodinger  equation.  We  shall  concentrate 
our  attention  mainly  on  this  computational  method. 

In  order  to  find  the  probabilistic  meaning  of  (l) , let  us  consider  the 
problem  of  a single  particle  performing  a random  walk  in  one  dimension  sub- 
ject to  the  following  conditions: 

A)  At  equispaced  time  intervals  TJ  the  particle  takes  a step  of 
length  h either  to  the  left  or  the  right,  with  equal  proba- 
bilities „ 

B)  If  the  particle  stays  at  a given  point  x during  a time  inter- 
val dt,  it  is  subject  to  destruction  there  with  probability 
V(x)dt  + o(dt) . (We  now  restrict  ourselves  to  the  case 

V(x)  * 0,  but  shall  later  drop  the  restriction.) 

If  we  let  hP(x,t)  represent  the  probability  that  the  particle  is  at 
the  point  x at  time  t (P(x,t)  is  thus  a probability  "density"),  clearly 
?(x,t)  must  satisfy  the  difference  equation 
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P(x,t)  - |P(x  - h,  t -X)(l  - 7(x  - h)X  + o(C)) 

+ fP(x  + h,  t - X)(l  - v(x  + h)t  + o(-C)) 

(We  have  merely  decomposed  the  event  of  being  at  the  point  x at  time  t 
into  the  two  alternative  ways  in  which  it  could  have  occurred*)  Subtract- 
ing P(x,t  - X)  from  both  sides  and  dividing  throughout  by  X}  we  obtain 

P(xJt)-P(x,t-'C)  _ 1 P ( x-h . t-X) -2P (x . t~X)+P ( h . t-ta 
X 2 % — “ 

- | P(x  - h,  t - X)v(x  - h)  - i P(x  + h,  t - r)7(x  + h)  + SEEL  . 

<-  1? 

The  left  side,  of  course,  is  a finite  difference  approximation  to  If 

we  now  require  that  the  space  and  time  increments  in  the  random  walk  be  re- 
lated by  the  equation  h2  ■ t;  f then  the  first  term  on  the  right  is  a second 

32P 

difference  approximation  to  * . Then,  in  the  limit  as  h 0 (and  hence 

^ yi 

^ °) 9 the  difference  equation  approaches  the  differential  equation  (l) . 

Thus,  equation  (l)  is  the  limiting  form  of  the  difference  equation 
describing  a random  walk  with  destruction,  where  P(x,t)dx  is  the  probability 
that  the  particle  is  between  x and  x + dx  at  time  t.  Equivalently,  the 
equation  describes  a diffusion  process  with  destruction,  where  P(x,t)  is 
the  particle  density  at  point  x at  time  t,  and  V(x)P(x,t)dx  is  the  time 
rate  of  destruction  of  particles  in  the  interval  (x  x + dx) . 

Under  the  former  interpretation,  if  we  start  the  particle  from  the 
origin,  we  must  require  that  our  solution  of  equation  (l)  satisfy  the  sub- 
sidiary condition 


■ 


. 


I . • • 


■ I *'  I I I I I I 1 1 1 

. X J 


. ' 

• • • 


, , • •••.  - 2 

- r r;  ' ' • ' ■ • ' ' 


H I I H III 


. 


. 


- h - 


(la) 


00 


for  all  t a o , j P(x,t)dx  * 1 


-oo 


and  the  initial  condition 


As  t 0+  , P(x,t)  -4  5(x)  . 


(lb) 


00 


(£>(x)  = 0 if  x / 0 . / 8 (x)dx  = 1) 

-oo 


If  V(x)  oo  as  x — > - oo  (we  shall  later  relax  this  restriction), 
the  solution  of  (l)  can  be  expanded  in  a series 

oo  “ X.t 

P(x,t)  = .2  c.e  3 ^-Cx) 

0=1  1 1 0 

where  the  X (>o)  and  are  the  eigenvalues  and  normalized  eigen- 

functions of 

2 

(2)  - Av|/(x)  = \ - v(x)  yyu) 

dx 


subject  to  the  conditions 

oo 


/ y (x)dxl  < 05  , / 


-oo 


oo 

> 

-00 


Y (*) 


'dx  < oo 


and  the  c.  are  so  chosen  as  to  satisfy  the  initial  condition  lb.  (If  we 
0 

look  for  a solution  of  (l)  in  the  fom  of  a superposition  of  functions  of 

-At  x 

the  type  ><pr(x)T(t)?  we  find  that  T(t)  * e and  that  the  space-depen- 

dent  part  vj/"  satisfies  equation  (2).)  The  eigenfunctions  ’vpr;(x)  form 

J 

an  ortho normal  set,  and  the  initial  condition  (lb)  yields  c.  = ^.(0). 

, 0 0 
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•We  note  here  that  (2)  is  Schrodinger1 s equation  for  stationary  states, 
where  A represents  E,  the  energy  of  the  state. 

Returning  to  the  limiting  form  of  the  random  walk  problem,  we  now  have 
the  result  that  the  probability  of  survival  of  the  particle  till  time  t is 

oo  -A  .t  oo  oo  - A .t 

(3)  / P(x,t)dx  - ? e 0 Yi(°)  f Aj/  (x)dx  s 2 b.e  0 

4o  0 4x>  J a*i  3 

If  we  arrange  the  Afs  in  order  of  increasing  magnitude,  the  first 
term  of  the  series  will  dominate  as  t — > oo  and  hence 


oo 


(a 


A = - lira 
1 t->oo 


log  f P(x,t)dx 


-00 


We  shall  employ  the  original  discrete  random  walk  to  obtain  an  esti- 
oo 

mate  of  J*  P(x,t)dx  for  large  finite  t.  Then  equation  (U)  enables  us  to 
-oo 


estimate  . 

Returning  to  the  discrete  random  walk,  we  inquire:  what  is  the  proba- 

bility p(x,9)  that  a particle  survives  the  time  interval  9(9  ^ X) , given 
that  the  particle  spends  the  entire  interval  at  the  point  x?  It  follows 
directly  from  our  description  of  the  random  walk  that  this  probability 
must  satisfy  the  differential  equation 


3p(*,9) 

Te 


V(x)p(x,9)  . 


Hence 


p(x,e) 


•v(x)e 


Now,  suppose  a particle  performs  the  discrete  random  walk,  starting 

from  the  origin  and  taking  steps  of  length  h * -7 as  at  time  intervals 

vn 
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X = 1.  Given  that  the  particle  follows  a particular  path  x(kt),  the  proba- 
bility that  it  survives  till  time  t(=nt^)  on  this  path  is 


nt-1 


(Ua) 


I I"  p(x(k*),X)  = exp  < - k|Q  V(x(k Z))Z 
k=o 


nt-1 


exp 


nt-I 

- 

k*o 


%(*(!£)  )ll  . 

n n | 

From  the  definition  of  the  random  walk,  it  follows  that  x(-)  * S,  v - , 

9 vn'  k n , 

where  is  the  sum  of  k random  variables  X-^Xg,  •••  each  X..  taking 
the  value  1 or  -1  with  equal  probabilities. 

We  have  computed  the  probability  of  survival  till  time  t on  a particu- 
lar path.  The  probability  of  survival  till  time  t,  irrespective  of  path* 
is  the  average  of  (Ua)  taken  over  all  2 possible  paths  x(kT?  . Since 
all  paths  are  equally  likely,  the  average  we  are  interested  in  is  the  simple 
arithmetic  average . We  symbolize  this  average  by 

(lib) 

Then,  for  large  n,  we 
(lie)  <^exp 

Thus,  by  equations  (U)  and  (l±c),  we  can  estimate  by 
(5)  X ~ -log  (ex?  {-I  At7^)) 


It  is  not  necessary  actually  to  compute  the  average  (Ub)  over  all 
possible  paths.  We  can  get  a statistical  estimate  of  the  average  by 
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considering  the  average  over  the  paths  m independent  random  walks. 

In  equation  (?)  we  have  now  derived  the  main  result  which  we  shall 
use  for  computational  purposes.  The  next  few  paragraphs  are  intended  to 
render  precise  some  of  our  remarks  in  later  sections,^ 

It  follows  from  the  central  limit  theorem  of  probability  that  in  the 

g 

limit  as  n 4 oo  and  k -4  oo,  the  distribution  of  Dk/\/n  becomes  Gaussian 

k 

with  mean  0 and  variance  - , 


i.e.  Prob 


{x<^  <x+  dx| 


V2lf  k/n 


exp 


2 
— x 

2k/n 


f dx 


Then,  as  n 4 oo  , the  distribution  function  of  the  random  variable 

i s r 

n k<nt  °k/v^i)  approaches  that  of  the  random  variable  J V(x(‘C))d‘£*  here 
the  path  x(T)  is  chosen  from  a space  in  which  the  measure  of  a set  of  paths 
is  the  probability  that  a particle  performing  a Gaussian  random  walk  will 
follow  a path  in  that  set,  (By  a Gaussian  random  walk  we  mean  a continuous 
random  walk,  with  independent  increments,  with  the  displacement  x(£)  at 
time  % governed  by  the  law 


(5a) 


Prob  y < x(T)  < y + dy  j-  - -/  - y2/2"Cdy  J.^>  . 

V(Sk//n)  <exp  - -J  V(x(C))dV 


Thus,  as  n -4  oo  , 


where  the  latter  average  represents  the  average  value  of 

exp  ^ - vfk  V(x(C“))dE^,  taken  over  the  space  of  continuous  paths  with  the 

^Tor  the  rigorous  development  of  the  theory,  see  [1], 
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measure  described  above.  This  average  is  just  the  Wiener  integral  of  the 
functional  exp  { - I v(x(t))drj.  , and  the  left  side  of  (5a)  msy  be  thought 
of  as  the  analogue  of  a Riemann  sum  approximating  the  Wiener  integral.  By 
(Uc)  and  (5a) 

t 


(6) 


| - J*  V(x('C))dt  Y>'C  P(x,t)dx 


For  the  purpose  of  computing  A^,  it  would  clearly  be  better  to  get  a 
sampling  estimate  of  the  right  side  of  equation  (5a)  than  of  the  left  side. 
However,  in  our  actual  sampling  process  we  cannot  deal  with  the  space  of 
continuous  paths  xtC) . We  can,  nevertheless,  speed  the  convergence  of  the 
distribution  of  - k<nt  ^(°k/\/n)  to  that  of  V(x(tO)(TE  by  taking  ^k/\/n 
directly  as  a Gaussian  variable  with  mean  0 and  variance  k/n,  rather  than 
waiting  for  k and  n to  become  large  enough  for  the  central  limit  theorem  to 
apply.  In  practice,  we  achieve  this  by  letting  5^  * X^  + * • • + X^,  where 
each  X^  is  Gaussian  (rather  than  discrete)  with  mean  0 and  variance  1. 

» Returning  to  equation  (5),  we  can  use  an  IBM  computer  to  find 
<\exp  f"n  klnt  v(Sk/\/n)}^  and  thus  X^,  as  follows: 

1)  Pick  values  for  n and  t (errors  involved  in  choosing  finite 
n and  t will  be  discussed  bel  ow) } 

2)  Choose  a random  number  from  a card  containing  Gaussian 
deviates  with  mean  0 and  variance  1; 

3)  Compute  V(-=)s 

Ml 

U)  Choose  a random  number  Xr>| 


. 

IM 
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Choose  a random  number  X, 
X-+...+X, 

Compute  V( p -)j 


,/5 


Continue  as  long  as  k < nt; 
Compute 


as  long 

i{v& 

n-l  - 


) . vffiS)  ♦ 


/ 


Compute  e 

We  have  now  found  one  term  in  the  average  ^ exp 
namely  the  term  corresponding  to  the  path  determined  by  our  particular 


, 

Vn  J 

{-  5 At  «7n’  } > ' 


choice  of  the  X^'s.  Each  choice  of  a sequence  X^,Xg,  •••  , corresponds 

-A 

to  a random  walk$  and  if  we  take  many  random  walks  and  average  e for  all 
the  walks , we  should  converge  to  an  accurate  value  of 

« kint  v^}>  • 

As  an  initial  test  of  the  method,  the  potentials  V(x)  = |x|  and 
2 

V(x)  * x (harmonic  oscillator)  were  studied  by  means  of  100  random  walks 
with  n = 100,  t^  = 3.7?,  tg  = ?.^  The  experimental  values  for  were  0.81 
and  0.7?,  compared  with  theoretical  values  0.81  and  0.71  respectively. 


Ill . ESTIMATION  OF  THE  ERROR  INVOLVED  IN  THE  METHOD 


The  feasibility  of  the  method  for  the  calculation  of  eigenvalues  of 
complicated  potentials  depends,  of  course,  on  the  accuracy  which  can  be 

*h?o  minimize  a numerical  error  discussed  in  Section  III  it  proves 
desirable  to  consider  two  different  va.ass  of  t.  For  a detailed  account 
of  the  experiments  on  these  potentials,  see  [2]. 
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obtained  in  a reasonable  number  of  random  walks.  There  seem  to  be  three 
major  possible  sources  of  error: 


A) 


B) 

0 


Error  involved  in  discretizing  the  problem  and  sampling  the 

1 °k\ 

values  of  the  finite  sum  - , 2 . V(-“)  instead  of  J V(x(r))dXr  # 

n n*t>  o 

i.e.  ^expj-i  <^exp{-f  V(x  (D)dT . 

Error  due  to  using  equation  (U)  for  finite  t. 


Sampling  error,  due  to  the  fact  that  we  sample  only  m paths 
out  of  the  space  of  all  possible  paths . 


Errors  B and  C,  which  we  shall  consider  first,  are  related.  It  will 
be  seen  in  the  consideration  of  statistical  errors  that  it  is  impractical 
to  make  t very  large.  Then  the  numerical  error  involved  in  using  equation 
(U)  for  finite  t is  twofold: 

log  b1 

l)  We  neglect  the  term  — ^ (see  eqn.(3)).  This  error,  which  may 

be  serious  for  moderate  t , is  easily  removed  by  considering  two 
different  times  t^  and  Then  it  follows  trivially  from  equa- 

tions (3)  and  (6)  that 

^ exp  ^ - I '1  V(x(Z))dZ 


(7) 


A - - 

i to-t 


log 


1 exp  j-  | 2 v(x(t-))dtr  ^ 

and  that  the  influence  of  b^  in  the  ratio  becomes  negligible  even 

for  moderate  t. 

In  the  expansion  of  J P(x3t)dx  we  omit  the  higher  exponentials 


in  t (see  eqn.(3)).  Ignoring  the  ration  l/b9,  we  see  that  the 


2) 


q n 
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. 
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-( A2-A  )t 

fractional  error  caused  by  this  omission  is  approximately  e ± . 

For  the  harmonic  oscillator,  A..  = , and  hence  the  fractional  error 

in  estimating  the  series  by  its  first  term  is  about  e For  t = b} 

this  is  less  than  0.5°/6  . 

The  sampling  error  C arises  from  the  fact  that  we  try  to  approximate 

the  exact  average  sT-u-  £ V(x(£)  )d£  (henceforth  called  > ) 

1 ® ^ *s 

by  the  average  - exp  -j  - ^ V(x^(tO  )dtjtaken  over  only  m independent 

paths  x^CC),  •••  y xm(tr) . The  expected  square  of  the  deviation  from 

<5**>  of  a single  value  exp  {-  l 7(xk(t'))dr}  is  given  by  < c*-2  > - <PC>^‘ , 

the  usual  expression  for  the  mean  square  deviation  of  a single  observation 

from  the  mean®  Then  the  mean  square  deviation  from  <oC>  of  the  average  of 

12  2 

m observations  is  - ( > - Thus,  the  relative  fluctuation  of 

m 9 

the  average  of  m observations  is  given  by 


i_  y«^>- < > 

Vfi  > 


1_ 


<c*>2 


- 1 


Fortunately,  we  are  able  to  estimate  < <*->*  for 


<oc2  > - <^expj^-2  J ?(x(t))dtj^)  . 

Then,  applying  the  results  of  equations  (3)  and  (6)  to  the  potential  2V, 


*LA.t  this  point  we  are  ignoring  the  discretization  error  A.  Other- 

n t 

wise,  we  should  properly  replace  J V(x^(V))dV  by  a Riemann  sum. 

o 


. 


' • . 

' T: 

Vl  *\  ...  -r  f li,  J&tlLtfJ 
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we  have 

^ .2.  . "Pf 

<.&**  > cL e 

where  is  the  lowest  eigenvalue  of  the  potential  2V;  ^-i(x)  is  the  lowest 

_ A™ 

eigenfunction  of  the  potential  2V;  « y^(°)  J <^(x)  dx. 

-oo 

The  relative  fluctuation  in  the  mean  is  thus 

__T 


(8) 


, 2 -2  A't 
bl  e 1 


- 1 


1 (2A1-1i1)t/2 

which,  for  large  t,  is  of  the  order  -7=,  e . For  any  positive 

vra 

1 

potential  V,  2 ^ 0,  and  hence  t must  be  kept  relatively  small, 

even  at  the  price  of  an  increase  in  error  B. 

In  the  case  V(x)  * x , we  know  a simple  change  of  scale 

yields  » 1.  Then  the  relative  fluctuation  in  <^J>  when  m - 100,  t = 5 
is  approximately 

i e(^-l)  5/2  „ 3o<£  . . 


This  result  is  not  so  disastrous  as  it  looks,  since  the  method  of  calcula- 
tion involves  the  logarithm  of  <©c>  and  is  not  so  sensitive  to  errors  in 

^“This  fact  is  easily  proved  by  the  variational  method.  On  the  other 

hand,  one  can  convert  our  reasoning  into  a probabilistic  "proof”  that 

2 A,  - \L  ^ 0 as  follows : if  u,  > 2 , then  for  sufficiently  large  t we 

-fLt  -2  2 2 

should  have  e x < < e and  thus  < oC  ^ — <<*■>  <03  but  this  would 

imply  that  for  sufficiently  large  t,  the  mean  square  deviation  of  a single 

observation  of  c*  from  the  mean  is  negative.  This  is  impossible. 
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<©^>  * Equation  (5)  says  A^  **  r log  <£p<>  . Then 


(9) 


dX 

/V 

X 

— d (log  <«*•>) 


1 

d<*> 

/■V 

tx 

<*<•> 

For  t = 5,  the  expected  percentage  error  in  our  determination  of  <oc> 
for  the  harmonic  oscillator  is  thus 


dX 

A 


iOI7  30C£  ~ 9%  • 


The  observed  relative  error  in  for  this  case  was  about  n- 

The  error  A is  not  easy  to  estimate.  What  we  should  like  to  find  is 
the  relative  error 


<^exp  j 

V(x('C))d'c|’^>  - <^exp  | 

; 5 klnt v(Sk/^] 

i> 

<^exp  -j 

f-  | V(xCC))dt 

}> 

W5 


A method  for  estimating  KTr  is  proposed  here,  but  not  rigorously  justified* 

In  the  case  of  the  harmonic  oscillator,  however,  we  are  able  to  calculate 

K 2 exactly  (Appendix  A) . The  answer  obtained  for  K 2 agrees  with  that 
x^  x^ 

predicted  by  our  more  general  method.  We  shall  discuss  the  method  briefly] 

it  is  worth  noting  that  in  practice  the  error  seems  unimportant. 

With  a given  continuous  path  x(^)  x^e  can  associate  another  path 

x^Ct),  which  is  a step-function  approximation  to  x("C) . We  define 


Then  ^exp 


{-* 


X (t)  - X (-) 

n n 


k<nt  ^ 


when  - ^ T*  < 

n n 


* <^exp  | 7(xn(t.))dt  ► 'y  . 


. 


■ 
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; ‘ - ' 
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,1 


By  definition  x^O  = x(-  [ntj,  and  thus  both  x(t)  and  x (z)  are  Gaussian 

random  variables.  If  we  plot  the  time  X versus  the  difference  ( X - -[nt^) 

n 

of  the  time  arguments  of  the  two  random  variables,  we  get  a sawtooth  func- 
tion which  grows  linearly  from  0 to  l./n  with  period  l/n.  On  the  average, 
the  step  function  lags  behind  the  continuous  function  in  time  by  an  amount 
l/2n.  Thus,  with  respect  to  functionals  which  depend  on  the  behavior  of 
x (r)  over  periods  of  time  large  compared  to  l/n,  we  may  replace  xn(£)  hy 
x(t*-  . (This  might  be  invalid  if  V is  very  rapidly  varying,  but  seems 

permissible  for  V(x)  * x . In  fact,  it  is  easily  verified  that  this  scheme 
gives  xnk('C)dr^> correctly,  to  the  first  order  in  l/n,  for  any  k.) 

We  therefore  conjecture  that,  to  the  first  order  in  l/n, 

f ^ • f 

"t  "fc 

<^exp  < - f V(xn(t))dt  - <^exp  < - ^ V(x(t- 


But  the  latter  quantity  is  the  same  as 
By  equations  (3)  and  (6),  for  large  t 


<r  { - ? r" 


V(x(t))dX  J^>  . 


t-  i. 

- f 2nV(x(t))dt 

O 


e 


- Ai(t-  k> 


Thus 


Kv- 


V -xi(t- 

— e 

=7£E 


^1 

- 

- 2n 


to  the  first  order  in  l/n«  For  n =*  100,  this  error  is  negligible. 


■[7]  * the  greatest  integer  less  than  or  equal  to  y. 
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Our  estimate  of  the  error  K^.  is  considerably  less  than  that  of 


Fortet  [3],  who  proves  rigorously  for  the  potential  V(x)  = x , that 

K 


2 t^j: 

2 4 5 V — 


X ' - ^ 

and  regards  the  right  side  as  a good  estimate  of  the  left,  Fortet !s  esti- 
mate follows  directly  from  his  preliminary  proof  that 


x2(X)dt^>  — <4: 


exp 


- f xn2(c)dr}  >1  ^ | 


£/2 

Vn 


This  result  is  true,  but  the  chain  of  the  inequalities  used  in  proving  it 
involves  some  very  wasteful  inequalities,  which  operate  to  make  the  right 
side  a very  high  bound  on  the  error.  The  highness  of  the  bound  is  most 
easily  seen  by  considering  the  case  of  fixed  n and  large  t.  Then  the  right 
side  approaches  + oo  as  t — > oo,  while  clearly  the  left  side  approaches  0 
(since  each  of  the  terms  on  the  left  goes  to  0) . 

We  prove  in  Appendix  A that  K — , correct  to  the  first  order 

1 ‘ x ^n 

in  (This  result  agrees  with  our  conjecture).  For  the  interesting  case 

of  n - 100,  t = 5,  we  thus  have 


K ' 


If  we  take  Fortet* s bound  on  K 2 as  an  estimate,  we  obtain  for  t 
n = 100 


- 5, 


K 2 - 16  - n.600^5 

Thus  the  latter  estimate  appears  too  high  by  a factor  of  ?000  and,  if 
correct,  would  indicate  that  the  method  is  totally  useless. 
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IV.  EXTENSION  TO  THE  HYDROGEN  ATOM 


The  harmonic  oscillator , of  course,  represents  one  of  the  simplest 
potentials!  the  results  obtained  in  the  study  of  the  hydrogen  potential 
are  of  more  significance  in  determining  the  real  value  of  the  method.  In 
Cartesian  coordinates  and  in  the  proper  units , Schrodinger!s  equation  for 


the  hydrogen  atom  is 

i2-ur  -2.,.  ^2 


1 f9~Y  ry  ry 


V X 


'of 


't  z 


Y - - > Y • 


x +y  +z 

We  encounter  here  the  following  complicating  features,  which  were  not 

present  in  the  case  of  the  harmonic  oscillator: 

1 1 


X.  V(r)  = - 


depends  on  three  coordinates.  Hence 


x +y  +z 

the  random  walks  must  be  performed  in  three  dimensions  rather 
than  one. 

2e  V is  negative,  has  a singularity  at  the  origin,  and  does  not  ap- 
proach oo  as  r oo  (If  we  want  to  continue  thinking  in  terms  of 
random  walks,  we  should  now  think  of  V as  representing  a proba- 
bility of  multiplication,  rather  than  destruction) * 

3.  The  spectrum  of  eigenvalues  is  partly  discrete  and  partly  con- 
tinuous, the  continuous  spectrum  arising  from  the  fact  that  V 
does  not  approach  oo  as  r *4  oo  9 The  eigenvalues  are  now 

\ * - rj  (n  * 1,  2,  3,  •••  ) and  the  entire  positive  real  axis, 

n 

Results  analogous  to  those  of  the  one-dimensional  case  can  still  be 
obtained.  Corresponding  to  equation  (3)  we  have 


. 
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l&M 


d-oo 
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p dtr 

J W)  \oo-X.t  oo  oo  \ , 

e 3 ^(0,0,0)  fff'jf 3(x,y,z)dx  dy  dz  + f e-A’tC(A)dA 


where  the  integral  is  the  contribution  of  the  continuous  spectrum.  For  con- 
venience, we  redefine  X as  i i and  can  now  write 

n <L  d 


where 


Cl  - Y^OjO.o)  jff  y (x,7,z)  dx  dy 

-00 


dz 


The  replacement  of  J ¥(r(t?))d'C'  by  a finite  sum  proceeds  as  before, 


72“’72T2” 

s,  +s,  s.  s, 

with  J — — — substituted  for  — . Here  S.  is  the  sum  of  k inde- 

n Vn  ^ 

pendent  and  identically  distributed  random  variables  X^,  •••  , X^.,  each  with 
mean  zero  and  variance  one  (similarly  and  S^)  • Thus  the  three-dimen- 
sional path  is  determined  by  three  independent  one-dimensional  random  walks. 

-t  A dt  . 


Then  - J*  V(r{X))&V  = J r(t7  by 


1 2 

n k<nt 


2“ — T"  mT~ 

sf  +sf  +sf 

kx  ky  kz 
n 


~ 2 


/ 


2 2 2 

sr+sf  +sf 

kx  ky  kz 


As  before,  we  approximate  X^  by  the  relation 


■-  . b/i  | 


, 
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</^p> 

<t^> 


Again  we  use  two  different  values  of  t to  eliminate  the  error  B, 

In  practice,  some  sampling  time  is  saved  by  the  following  scheme.  We 

note  that  x(ZT) , being  Gaussian  with  mean  0 and  variance  IT,  has  the  same 

distribution  as  \/t  x(£) , Therefore  r(t»)  has  the  same  distribution  as 

t t 

Vt  r(^) . It  follows  that  J — d'T  has  the  same  distribution  as 
1 i i r J 

vt  J'  rfgy  ^ • Then,  what  we  actually  sample  is  the  random  variable 


c<- 


n 1 

-J  ?5tT 


dtr 


100 


/ 


— 2 — T 

sf  +sf  +sf 

kx  ky  kz 


We  then  compute 


\ by 


the  relation 


t2“tl 


log 


</^> 


The  values  of  t used  must  be  large  enough  to  make  error  B (discussed 
in  section  HI)  small*  The  magnitude  of  this  error,  as  we  have  already 
seen,  is  about  e (the  sign  of  the  exponent  differs  from  that  in 

the  similar  expression  for  the  harmonic  oscillator  because  the  solutions 
now  increase,  rather  than  decrease,  exponentially).  For  the  hydrogen 
atom,  - A 2 * 0,2?.  The  corresponding  quantity  for  the  harmonic 
oscillator  is  l.U.  When  t^  * 9,  the  fractional  error  in 


1 ■ 
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\e  ) due  to  B is  thus  about  10%  , and  the  consequent  fraction- 
al error  in  about  2°/0 . For  t2  23  16,  the  accuracy  is  better. 

What  happens  to  the  sampling  error  C when  we  use  such  large  values  of 
t?  If  we  sample  m paths,  we  recall  that  the  relative  fluctuation  of  the 


X m o rk(*)dt'  1 (fI1-2*1)t/2 

average  - 2 e is  approximately  -7=  e $ u-,.  the  lowest 

m k=±  vm  1 

2 

eigenvalue  of  the  potential  - - , is  easily  calculated  as  2.  Thus 

Pi  - 2^i  - 1,  The  corresponding  term  in  the  statistical  error  computation 

for  the  harmonic  oscillator  was  .I|lU,  Hence  the  statistical  fluctuation 

for  a given  number  of  samples  m and  given  t is  considerably  larger  for  the 

1 2 

potential  - - than  for  the  potential  x „ 

Intuitively,  one  can  explain  this  increase  in  the  statistical  error  in 
the  following  manner:  The  eigenvalue  of  the  hydrogen  potential  depends 

mainly  on  the  behavior  of  the  potential  in  the  small  region  near  the  origin. 
We  see  this  fact  immediately  when  we  note  that  the  main  contribution  to  the 


<• 


J* 


t 1 


ritry 


dr 


average  V e >>  comes  from  the  few  paths  r Cc)  which  stay  close  to 

the  origin,  Plence,  in  order  to  calculate  the  eigenvalue  accurately  we  must 
get  a detailed  account  of  the  behavior  of  the  potential  in  the  critical 
small  region  about  the  origin.  Naturally,  many  samples  are  required  before 
we  get  a sufficient  number  of  paths  staying  close  to  the  origin  to  sample 
that  region  accurately.  Equivalently,  if  we  do  not  take  a very  large  num- 
ber of  samples,  there  is  a good  chance  that  the  few  critical  paths  which 
stay  near  the  origin  will  not  turn  up  In  the  proper  frequency,  thus 
causing  a large  error  in  the  eigenvalue , In  the  case  of  the  harmonic 
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oscillator,  however,  there  is  no  such  small  critical  region  which  must  be 
carefully  sampled*  Remarks  similar  to  these  about  hydrogen  may  also  be 
made  about  atomic  potentials  in  general . 

Returning  to  our  computations,  we  see  that  for  900  samples  with  t = 9 


m 


the  expected  relative  fluctuation  of  - klie 


•t  i 

i ^ 


dt 


V 900 


9/2  ~ 35  ~ 300°^,  . 


is  about 


If  we  compute  more  carefully,  using  the  exact  expression  for  the  relative 
fluctuation  given  in  equation  (8),  we  obtain  100°^  instead  of  300^.  It 
follows  that  for  m = 100  and  t = 9, 


d X 

x 


1 d <£*-> 

t \ 


1 

TZ5 


100°7o 


22®?,  . 


For  t = 16,  the  error  is  of  course  much  greater.  For  the  initial  group  of 
900  hydrogen  random  walks,  with  t^  = 9 and  t^  = 16,  computations  according 
to  equation  (10)  yielded 


.503  (true  value  *p00) 


In  view  of  the  preceding  remarks,  this  accuracy  was  quite  surprising. 
An  examination  of  the  data,  however,  seems  to  confirm  our  previous  analy- 
sis. Of  the  900  oaths  sampled,  I4.  yielded  extremely  large  values  of 

-1  x 

J rTFT  ^ * without  these  I4.  paths,  which  clearly  cannot  be  counted  on  to 
recur  regularly,  the  average  e4 

the  computed  value  for  A about  20°/o  smaller  (note  that  the  logarithmic 
method  of  calculation  is  not  too  sensitive  even  to  a large  fluctuation  in 

6Uoc>). 


would  have  been  smaller,  and 
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Thus,  there  seem  to  be  good  grounds  for  attributing  the  high  accuracy 
in  this  determination  of  A^  to  luck.  Nevertheless,  it  seemed  worthwhile 
to  run  another  group  of  samples.  This  was  easily  done  in  connection  with 
experiments  on  helium  (described  below),  since  each  helium  random  walk 
yields  two  hydrogen  random  walks  as  a by-product.  Two  more  groups  of  900 
hydrogen  walks  were  thus  obtained. 

The  qualitative  appearance  of  the  data  was  similar  to  that  of  the  first 

group  of  walks.  In  the  second  group  of  900,  more  than  Q0°/o  of  the  contri- 
/ L°£,V 

bution  to  C e y came  from  three  walks  yielding  high  values  of  ot  , and 
the  calculated  value  of  was  0.36.  If  a few  more  high  values  of  had 
been  recorded,  the  calculated  value  of  A^  would  have  been  much  closer  to 
0.50j  but  this  again  would  have  been  purely  a matter  of  luck.  In  the  third 
group,  there  were  only  two  large  values  of  o(,  with  the  result  that  the  low 
value  * 0,29  was  obtained. 

A more  detailed  analysis  of  some  aspects  of  the  hydrogen  random  walks 
is  given  in  section  V and  in  Appendix  3. 

The  preceding  discussion  of  hydrogen,  which  may  be  extended  to  other 
atomic  potentials,  makes  one  fact  rather  clear;  for  t great  enough  to  make 

our  numerical  approximations  good,  the  statistical  fluctuation  in 

dt 

/ £ Vpf  \ 

(e  K / becomes  prohibitively  large  unless  we  use  a huge  number  of 
samples.  Equivalently,  our  computed  value  for  is  determined  almost  en- 
tirely by  a very  small  group  out  of  the  900  random  walks.  Hence,  most  of 

I 

the  data  is  wasted  and  our  large  group  of  samples  acts  statistically  like 
a small  group. 
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There  now  arises  the  ^question  of  whether  there  is  any  other  method  of 
utilizing  our  data  to  find  A_^,  by  means  of  some  computation  less  subject 
to  statistical  fluctuations.  In  the  next  section  we  shall  deal  with  such 

t 

a method,  which  we  call  the  slope  method,  as  distinguished  from  our  previous 
mean  value  method. 


V.  THE  SLOPE  METHOD 5 RESULTS  FOR 
HYDROGEN,  HELIUM,  AND  n DOUBLE  HYDROGEN’1 

The  method  for  finding  A ^ which  is  to  be  presented  here  has  the  fol- 
lowing advantages  over  the  mean  value  method: 

1)  It  makes  use  of  the  finer  structure  of  the  distribution  of  the 

;t  d 

IrTiET  anc^ 

o ^ ' 

2)  It  is  quite  insensitive  to  erratic  individual  observations  of 
the  random  variable*  — - 

Consider  the  random  variable  r^TfJ  * this  random 

variable  assumes  depends  on  the  particular  path  r(r)  which  is  chosen  for 
the  path  of  integration.  We  have  already  observed  (section  IV)  that 
J1  ^as  the  same  density  function  as  ]/i  J* 


«1 

Thus,  if  we  let  /?(x^  be  the  density  function  of  j 

r i -t  0 

(i.e.  />(x)dx  = Prob  -I  x < J*  ^ 

L o 

t J i 

p dr 

<J  , , 


Hr) 


< x + dx  f ) we  have 
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P dr 

\ft  J vi'Z  ) 


J'  X /(x)  dx  ~ 
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A-,t 


bn  e 


t—^oo  1 


■ 


' 


■ 


« 


/ 


. 

c ti  b ■ n f ' .fSJi  l 

' 

' ■ • 


. ' 


■ 


- 23  - 


or,  writing  t for  t,  we  have 

oo 
o 


(n) 


J etX/’(x)d x 


t 00 


V 


Alt' 


One  can  verify,  by  a simple  integration,  that  this  asymptotic  behavior  of 


oo 


tx 


J e y?(x)dx  can  be  achieved  by  setting 


(12) 


/(x) 


VhX1 


X — 00 


2 \TWT. 


By  an  unproved  Tauberian  theorem  which,  however,  is  intuitively  reasonable 
and  similar  in  form  to  the  Hardy-Littlewood  Tauberian  Theorem,  this  is  the 
only  asymptotic  behavior  of  ^(x)  which  will  make  equation  (ll)  true.  A 
result  similar  to  this  can  be  obtained  for  any  homogeneous  potential  func- 
tion V. 


Now  we  can  determine  A as  follows: 

1 1 


dZ 


l)  Find  the  density  function  of  J*  by  replacing  the  integral 


with  a finite  sum  and  performing  random  walks  as  before,  to 

get  sample  values  of  P • 

b ® 2 

2)  Letting  x = J1  \ , plot  log  /’(x)  vs.  x for  large  x.  It  is 
o ^ 

to  be  hoped  that  a fairly  straight  line  will  result.  The  slope 
of  this  line  should  be  - ♦ 

First  we  present  uncritically  an  example  of  how  this  technique  is  ap- 
plied to  the  hydrogen  potential. 


’The  theorem  can  be  proved  in  idle  integral: |^form 
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Clearly,  we  cannot  experimentally  find  a continuous  density  function 

1 dZ 

/?(x)  . Our  experimental  data  consists  of  900  values  of  J1  obtained 

o rv  / 

from  900  random  walks.  The  closest  we  can  come  to  plotting  a density  func- 
tion is  to  plot  a histogram,  which  shows  the  number  of  occurrences  of 

1 dZ 

S rTcT  each  of  a sequence  of  intervals  of  length  A;  e.g.,  if 

o ' ' 

we  take  A = 0.1,  we  obtain  the  histogram  of  Graph  1.  The  bar  covering  the 

/I  j 

'-T'T 

o ' ' 

lay  between  1.5  and  1.6  (=  1.5  + A).  Of  course  the  choice  of  A does  much  to 
determine  the  appearance  of  the  histogram;  a small  A gives  an  irregular 
curve,  while  a large  A obscures  the  structure  of  the  histogram  by  "smooth- 
ing” it  too  much.  Since  the  purpose  of  our  histogram  is  only  to  study  the 
behavior  of  the  tail,  it  is  clear  that  A = 0.1  makes  the  tail  too  irregular 
for  reliable  use  in  our  calculations.  If  we  replot  the  histogram  with 
A = 0.2  (graph  2),  the  tail  is  considerably  smoother  but  still  retains  some 
fine  structure. 

2 

Now,  if  we  plot  log  f{ x)  vs . x for  large  x,  we  expect  a straight 

line  with  slope  - • faking  values  of  f from  graph  2 and  plotting 

2 ^ 

them  against  x on  semi-log  paper,  we  obtain  graph  3*  Here  the  "tail”  of 
p is  considered  as  beginning  at  x * 1.2,  just  to  the  right  of  the  peak 
of  f . The  question  immediately  arises:  what  is  the  best  line  to  draw 
through  the  points?  We  reason  as  follows  (this  argument  and  the  general 
character  of  the  graph  are  typical  of  any  potential): 

1)  The  first  few  points  at  the  left  are  not  really  in  the  tall 
since  x is  not  large  enough.  This  is  also  clear  from  the 
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graph,  where  it  is  seen  that  the  points  cannot  be  well  fitted 

by  a line,  i.e.  they  do  not  lie  in  the  region  where 
— x2/4^2_ 

P(x)  ~ const,  e 

2)  The  larger  the  x-coordinate  " int,  the  more  valid  is  the  ap- 


Thus,  the  data  points  becpme  statistically  less  reliable  as  the  numer- 
ical approximation  improves.  Consequently,  we  reason  that  the  best  line  is 
the  one  going  through  a set  of  points  whose  x-values  are  as  large  as  pos- 
sible, but  which  still  represent  enough  data  to  be  fairly  reliable.  Of 
course,  the  test  of  the  "reliability”  of  a point  is  that  it  lie  on  a line 
with  several  other  points.  Near  the  mid-portion  of  the  graph  there  gener- 
ally appears  a sequence  of  such  points,  followed  by  a number  of  "wild" 
points  which  represent  insufficient  data.  Drawing  a line  on  graph  3 ac- 
cording to  these  principles  and  measuring  its  slope,  we  obtain  A = .1*8 . 
Obviously  it  is  difficult  to  give  a fair  estimate  of  the  error  in  such  a 
method.  A good  upper  bound,  however,  is  obtained  by  looking  at  the  slopes 
of  the  most  extreme  lines  one  might  attempt  to  draw  through  the  data 
points,  and  seeing  how  much  those  slopes  differ  from  that  of  the  original 
line.  In  this  case  the  answer  is  quite  certainly  less  than  1096. 


proximation  ^(x)  ^ const,  e . But 


3)  the  points  which  represent  really  large  values  of  x represent  very 
little  data,  i.e.  there  are  few  occurrences  of  very  large  values 
. Hence  the  points  with  the  largest  x-coordinates  are 
subject  to  large  statistical  errors. 
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"The  discussion  gains  "ibility  who..  we  consider  the  results  for 

the  hydrogen  potential  obtained  from  1800  mere  random  walks.  The  histogram 

(graph  U)  has  the  same  shape  as  graph  2 and  has  a peak  at  the  same  value 

of  x.  When  we  plot  log  />(x)  vs.  x for  large  x,  we  get  graph  5>,  which  shows 

a marked  improvement  over  graph  3.  With  the  exception  of  the  last  two  points 

on  the  right,  which  represent  insufficient  data,  the  points  quite  certainly 

lie  near  a line.  Calculating  the  slope  of  a line  fitted  by  eye,  we  find 

A = .52 . Here  the  error  seems  less  than  sh°. 

It  is  worth  noting  that  about  30°/o  of  the  random  walks  yielded  values 

-y  \ which  are  large  enough  to  be  considered  in  the  "tail11  of  the 
o ^ ' 

distribution.  Thus,  in  the  case  of  the  hydrogen  atom,  one  can  say  roughly 

that  the  slope  method  utilizes  30°7o  of  the  data  while  the  mean  value  method 

t 

makes  use  of  only  about  19&  (since  the  average  e°  ^ is  determined 

largely  by  a handful  of  walks  yielding  large  values  of  the  integral) . This 
difference  in  efficiency  of  the  two  methods  seems  critical  in  the  case  of 
hydrogen,  since  the  mean  value  method  is  infeasible  for  that  potential 
T<rhile  the  slope  method  is  reasonably  accurate.  We  can  get  an  idea  of  the 
limitations  of  the  slope  method  by  considering  a more  difficult  potential. 

In  the  proper  units,  the  stationary  SchrSdinger  equation  for  helium 
can  be  written  as 


| q2  'k  (5?,  + \ V k 
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In  order  to  find  X ^ by  the  slope  method,  we  must  determine,  by  sampling, 
the  distribution  of  the  random  variable 

/ f-J—  + —2 i— 

where  r^(£)  and  r^(t)  are  the  paths  of  two  independent  Gaussian  random 
walks.  As  before,  we  discretize  the  walks  and  sample 


.1 


100 

k*l 


2 


/ 


2 2 2 
S.  +sf  +sf 
kx2  ky^  kz2 


(1U) 


* 

► 


^ <VV*+(VV*+(VS)2  - 

,Sky1,SkZi 

dependent  three-dimensional  random  walks,  as  in  Section  IV. 

Before  we  discuss  the  results  for  helium,  let  us  consider  a slightly 
simpler  potential,  whose  lowest  eigenvalue  can  be  calculated  as  a by-pro- 
duct of  the  helium  computation.  Consideration  of  this  potential  will 
bring  out  the  difficulties  involved  in  the  multi-electron  potentials.  The 
potential  we  deal  with  is 


where  (S. 


kx1 


) and  (S^.  ,Skz  ) determine  the  paths  of  two  in- 
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which  corresponds  to  a helium  atom  without  interaction  between  the  elec- 
trons. This  potential  has  been  christened  "Double  Hydrogen".  The 
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Schrodinger  equation  for  double  hydrogen  is 


I «-»  — v \ . ± n 2 . j ft 


(15')  i V ‘Xf  (q,r2)H  | C’fy,r2)  + (• 
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j£_) 

>— ■*  1 ' 

ir^  i 


Ub  ,v9)  « - +\r- 


It  is  easy  to  show  that  the  eigenvalues  of  this  equation  are 
'*■ = Pi + jv  where  p.  and  p^  are  eigenvalues  (not  necessarily  distinct) 


or 


(16) 


' £ 

The  eigenvalues  of  (l6)  are  \i^_  * - (k  - 1,2,3,  •••  ) and  the  entire 

k 

positive  real  axis.  The  lowest  eigenvalue  of  (l5)  is  thus  = -Ii,  and 
the  second  lowest  is  Ap  = -2.5.  As  before,  we  now  redefine  A so  that 
, ^ = U,  “ 2.5.  (It  also  follows  simply  from,  our  probabilistic  ap- 
proach that  A1  « 2ytj , since 


t 

X 

o 1 


, /(nfw  + ^m)ax 


> 


X ^ lQg-S  e 

1 


but  r^(T)  an(^  t2 (^)  are  the  paths  of  independent  random  walks,  and  hence 

e e e . 

The  desired  result  follows.) 

Of  course,  since  we  know  that  A_  - 2ii_  , one  statistical  method  for 

1 : JL 

finding  would  be  to  find  by  saiy  rag  the  potential  - , and  then  t< 
double  the  result.  This  scheme,  however,  breaks  down  as  soon  as  the 


-4 


potential  contains  an  interaction  term  between  it.  and  r0«  Hence,  since 
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we  want  to  test  the  general  applicability  of  the  method  to  multi- electron 
potentials,  we  do  not  employ  this  short-cut,  but  sample  the  random 
variable 


Plotting  the  logarithm  of  the  density  function  of  this  random 

2 

variable  versus  x , we  obtain  graph  6.  The  graph  is  based  on  900  values 

1/2  2 \ 

of  S (r  Wj  + r ~'(v) ) ^ 9 °“  &bout  500  are  represented  on  the  graph 

(the  other  values  were  too  small  to  be  in  the  "tail”  of  /’(x)).  On  the 
basis  of  the  solid  line  drawn,  we  find  ^ = 5.6  (true  value  = U.0) . It 
is  evident  that  the  data  points  do  not  clearly  determine  a line.  If  we 
ignore  the  points  on  the  left,  arguing  that  they  are  not  yet  in  the  tail, 
and  draw  the  dotted  line  based  on  the  right  points,  we  get  A^  = U.6.  This 
result,  however,  is  hardly  reliable.  If  we  increase  A from  0.6  to  0.8,  we 
obtain  graph  7,  which  is  considerably  smoother.  The  slope  of  the  line 
through  the  first  four  points  on  the  left  yields  = 8.5.  Clearly 
these  points  are  not  yet  in  the  tail.  The  next  three  points  represent 
only  29*  8,  and  6 walks  respectively  and  are  subject  to  too  great  statis- 
tical fluctuations . The  expected  statistical  fluctuation  of  a point  rep- 
resenting n walks  is  about  \/n.  Hence  the  limits  of  inaccuracy  cf  the 
data  points  are  represented  approximately  by  the  vertical  bars  we  have 
drawn.  We  see  that  the  points  are  not  inconsistent  with  the  true  value 
A = U,  but  they  hardly  give  us  any  further  information. 
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The  redeeming  feature  of  the  slope  method,  in  the  case  of  the  hydro- 
gen  potential  (and  thus  also  for  the  potential  -)  was  that  the  exponential 
” tail"  of  the  distribution  started  early  enough  so  that  the  points  in  the 
tail  represented  enough  walks  to  be  statistically  reliable.  This  feature 
does  not  seem  to  hold  for  double  hydrogen.  We  now  proceed  to  a more  pre- 
cise discussion  of  this  point. 

A simple  extension  of  the  reasoning  leading  to  equation  (12 ) yields 


(17)  r(x)  ~ — — 
-WCD  2/WX 


-x2/Ut 


2\ZrirX 


a ”v " ^ 2 

Thus,  an  estimate  of  the  goodness  of  the  approximation  (12)  is  given  by 

% 

the  ratio  of  the  first  two  terms  on  the  right  side  of  (l7)j  this  ratio  is 
approximately 

2/1  1 \ 
x (ur  up 


When  x is  large  enough  to  make  this  ratio  large,  we  can  say  we  are  in  the 
tail®  Thus,  the  criterion  for  being  in  the  tail  is 

(18)  tj  (t2  ~ a~}  58  M 

where  M is  some  constant,  whose  value  depends  on  how  good  we  want  the 
approximation  (12 ) to  be. 

1 11 

For  the  hydrogen  potential  -,  we  have  -r-  “*  “ 8-2  = 6,  whereas  for 

r Ar)  A-i 

1111  <-  ± 

double  hydrogen  " j;  “ . Thus,  if  the  tail  for  hydrogen 

starts  at  xq,  the  tail  for  double  hydrogen  starts  at  \/IJ 0 Now  we  ask 
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the  question:  how  mar$r  random  walks  for  double  hydrogen  must  we  perform 

in  order  to  obtain  an  accuracy  comparable  to  that  obtained  for  hydrogen? 

If  we  perform  random  walks  for  hydrogen,  the  number  of  walks 

yielding  values  of  f in  the  tail  of  /H(x)  is  about  J /^(x) dx. 

n x 

o 

To  achieve  a comparable  accuracy  for  double  hydrogen,  we  want  to  get  the 
same  number  of  points  in  the  tail.  Hence  we  require 

00  pCD 

/ fn{x)  & = "dh  J^r  ^dh(x) 


n 


X 


where  = number  of  double  hydrogen  random  walks 


/?I)H(x)  = density  function  of  J* 


+ wfrj  dV  • 


Ignoring  the  constant  factor  in  equation  (12),  we  thus  have  the  order- 


cf -magnitude  equation 
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Using  the  asymptotic  series  for  the  two  integrals,  we  have,  if  xq  is  large. 
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Looking  at  graph  3,  we  see  that  the  tail  seems  already  to  have  started  at 

2 

x =2.  The  above  calculation  is  inaccurate  for  such  a small  value  of 
o 

(we  could  calculate  ^DH/N^  more  exactly,  but  it  hardly  seems  worthwhile); 
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it  yields  the  qualitatively  correct  result  that,  for  a given  accuracy  in 


the  lowest  eigenvalue,  we  require  many  more  samples  for  double  hydrogen 
than  we  do  for  hydrogen. 

Formally , the  reason  for  the  increase  in  the  number  of  samples  re- 
quired is  the  creeping  together  of  the  reciprocals  of  the  first  and  second 


gen,  we  sample  a single  random  variable  X*  for  double  hydrogen  we  sample 

the  sum  X + Y,  where  X and  Y are  independent  and  identically  distributed. 

Roughly  speaking,  if  p is  the  probability  that  a value  of  X lie  in  the 

tail  of  the  distribution  of  X,  then  the  probability  that  a value  of  X + Y 

2 

lie  m the  tail  of  the  distribution  of  X + Y is  only  about  p . This  argu- 
ment can  be  made  quantitative  and  leads  to  our  previous  result. 

We  conclude  that  900  random  walks  constitute  far  too  small  a sample 
for  the  double  hydrogen  potential  or,  more  generally,  for  ary  multi-electron 
potential.  The  900  helium  random  walks  confirm  this  conclusion.  The  loga- 
rithmic histograms  with  A = O.U  and  A = 0.8  respectively  are  given  by 
graphs  8 and  9.  The  graphs  are  consistent  with  the  true  value  X - 2.90, 
but  not  enough  points  are  in  the  tail  to  allow  an  accurate  determination. 

In  the  next  section  we  discuss  some  schemes  for  making  our  method 
feasible  for  the  study  of  more  complicated  potentials.  One  more  point  is 
worth  mentioning,  however,  in  connection  with  double  hydrogen.  We  plotted 
a logarithmic  histogram  (graph  7)  for  this  potential,  based  on  900  values 


eigenvalues  (see  eqn.  (18)).  n Physically” , the  reason  is  this:  for  hydro- 


Let us  denote  these 


900  values  by  An  +-  A2  (1=1,  o«*,900)  where 
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However,  since 


/ 


o 


1 


are  independent,  A^  + A, 

i 


2 (i  ^ j)  or 
a 


A-,  + A-.  (i  / j)  or  A0  + A*  (i  / 3)  also  constitute  legitimate  sample 

JL  ■>  1 » tma  + €a.  % 
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values  of  the  random  variable  X.  Thus,  out  of  our  1800  values  A^  , , 

i i 

we  can  manufacture  1800  x 1799/2  **  160,000  values  of  the  random  variable  X, 
instead  of  merely  900  values*  Of  course, these  160,000  values  are  not  all 


a maximum  of  only  900  independent  values  of  X in  the  set  of  160,000.  How- 
ever, one  feels  intuitively  that  the  dependence  is  not  too  strong  and 
that  something  might  be  gained  by  using  the  additional  data.  (The  whole 
situation  can  be  treated  quantitatively,  but  we  omit  the  treatment  here. 

An  analogous  problem  is  that  of  sampling  the  distribution  of  the  sums 
rolled  on  a pair  of  dice.  We  roll  900  pairs,  but  then  combine  the  1800 
single  rolls  into  all  possible  pairs.) 

Actually,  we  did  not  combine  the  data  into  all  160,000  possible  pairs, 
but  considered  only  18,000  pairs.  Graph  10  is  the  logarithmic  histogram 
based  on  these  pairs.  The  improvement  over  graphs  6 and  7 is  marked,  and 
we  claim  about  6°/o  accuracy  for  the  result  X * 3.99.  One  might  tend  to 
challenge  this  result  on  the  ground  that  we  are  getting  n something  for 
nothing”,  and  thus  violating  what  some  consider  to  be  a fundamental  law  of 
nature.  However,  we  should  remember  that  the  data  already  contained 
enough  information  to  find  the  lowest  eigenvalue  of  ” single”  hydrogen  to 
about  5%  (graph  5),  and  thus  (by  our  remarks  on  p,  28)  the  lowest  eigenr- 
value  of  double  hydrogen  to  the  same  accuracy.  Graphs  6 and  7 merely 
represented  inefficient  use  of  the  data. 
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This  device  of  "multiplying"  the  data  fails,  of  course,  when  an  inter- 
action terra  is  present  (as  in  the  case  of  helium)  . Our  final  section  deals 
with  procedures  feasible  for  this  case, 

VI.  TECHNIQUES  FOR  SPEEDING  THE  SAMPLING  PROCESS ; 

HIGH-SPEED  COMPUTING  AND  IMPORTANCE  SAMPLING 

The  experiments  described  in  the  previous  sections  were  perfomed 
mainly  out  of  curiosity,  to  see  whether  equation  (5)  is  verifiable  within 
a reasonable  number  of  samples.  At  a later  stage,  one  is  naturally  led  to 
ask  whether  this  sampling  method  is  of  practical  use,  e.g.  as  a possible 
competitor  with  the  variational  method  for  finding  the  lowest  eigenvalue. 

The  authors  are  unaware  of  just  what  potentials  are  of  real  interest  but 
unamenable  to  variational  treatment.  However,  there  are  a number  of  tech- 
niques which  should  make  the  present  scheme  feasible  for  the  treatment, 
without  excessive  sampling  time,  of  any  system  whose  dimensionality  is 
small.  These  techniques,  which  will  be  dealt  with  briefly  here,  are 

1)  The  use  of  a high-speed  digital  computer,  eliminating 
entirely  the  "slow"  punched-card  operations; 

2)  The  employment  of  sophisticated  sampling  schemes,  which 
we  call  "importance  sampling" . 

Both  techniques  can  be  used  together,  effecting  a further  reduction  of 
sampling  time , 
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For  an  indication  of  the  times  involved  in  the  use  of  IBM  methods,  we 
note  that  900  hydrogen  random  walks  of  100  steps  each  require  about  90  IBM 
hours.  To  accumulate  the  same  data  on  the  National  Bureau  of  Standards 
Western  Automatic  Computer  (SWAC),  about  one  hour  is  required.  No  hydrogen 
walks  have  yet  been  performed  on  the  SWAC,  but  harmonic  oscillator  results 
will  soon  be  available.  One  can  say,  in  general,  that  for  our  purposes  the 
SWAC  is  faster  than  IBM  by  a factor  of  at  least  5>0 . The  programming  of 
these  random  walk  routines  is  quite  elementary,  involving  mainly  the  itera- 
tion of  simple  operations . The  problem  of  supplying  random  numbers , to 
govern  the  random  walks  is  of  some  interest.  External  sources,  such  as 
prepared  tables,  are  unfeasible,  since  the  machine  must  generate  its  own 
random  numbers  if  it  is  to  operate  at  electronic  speeds.  A number  of 
arithmetic  schemes  for  the  internal  generation  of  "random”  numbers  have 
been  proposed  (the  numbers  so  generated  are  only  pseudo -random,  since  they 
are  completely  determined  by  the  numbers  inserted  into  the  arithmetic 
scheme),  and  a particularly  simple  method  of  Dr.  J,  B,  Rosser  has  been 
found  to  produce  digits  which  satisfy  the  usual  tests  of  randomness. 

Regardless  of  whether  SWAC  or  IBM  equipment  is  being  used,  one  can 
gain  another  factor  of  10  or  more  by  the  use  of  importance  sampling.  This 
topic  has  been  discussed  in  recent  literature  (see  [5])  and  is  of  enough 
simplicity  and  interest  to  warrant  a brief  review. 

Suppose  wTe  have  a random  variable  X which  is  continuously  distributed 


b 


interval  (a,b)  with  the  probability 
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over  some 
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' and  suppose  we  are  interested  in  finding  the  expected  value  f of  some 

b 

function  f (X)  • Then  f = J f (x)/?  (x)dx.  (Alternatively,  and  for  greater 

a 

clarity,  we  sometimes  write  f = E^?(f(X)),  i.e.  f is  the  expectation  of 
f (X)  over  a population  in  which  X has  the  density  function  /?(x) .)  The  in- 
tegral can,  of  course,  be  estimated  by  sampling:  we  merely  pick  many  values 

of  X from  a population  governed  by  the  density  function  /?(x),  compute  f(X) 
for,  each  value  of  X,  and  average  the  results.  Suppose  we  also  compute 


(f  - f )d  ~ J*  (f(x)  “ f)2  /?(x)6x  = (Ty  p . 
a 

Then  it  is  well  known  that  the  expected  error  in  our  estimate  of  f , if  we 
take  as  our  estimate  the  average  of  n sample  values,  is  -i-  * • 

■s/rT  s>f 

Importance  sampling  is  based  merely  on  the  recognition  of  the  fact  that 

it  is  possible  to  bias  the  sampling  procedure  in  such  a way  as  to  leave  f 

unchanged  but  to  make  smaller.  Specifically,  suppose  y(x)  is  another 

probability  density  function  defined  on  the  interval  (a,b)  with 
b 

^(x)dx  * 1.  We  certainly  have 
a 

(19)  e p (f(x) ) = y3  f(x)/(x)dx=y  (x)dx  , 

Thus,  in  our  notation, 


£U)P( 
w (x) 


00 


E p (f  (X) ) = E ^ 

i.e.  the  expectation  of  f(X)  over  a population  governed,  by  the  density 
function  ^ is  the  same  as  the  expectation  of  f f/(f  over  a population 
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governed  by  the  density  function  The  interpretation  of  this  fact  is 


f ’ 


simple:  we  bias  the  density  function  of  the  population  by  the  factor 

P 

and  hence  must  ’’weight"  our  values  of  f by  the  factor  if  we  are  to  get 

the  correct  answer.  Let  us  compute  the  expected  error  in  the  average  of  n 

_ ]_ 

samples,  if  we  estimate  f by  this  scheme.  The  error  is  — „ . where 


oy, 

i, 


dx 


By  the  proper  choice  of  the  function  y(x),  ^ may  be  made  much 

smaller  than  - . In  particular,  if  we  knew  T (which  is  of  course  what 
we  are  trying  to  find),  we  could  let 


(20) 


(^(x)  = 


p{-x)  f(x) 


and  thus  have  no  error  in  our  estimate.  In  practice,  we  do  not  know  f,  and 
it  is  not  easy  to  concoct  a population  governed  by  an  arbitrary  density 
function  fU).  We  can,  however,  use  whatever  we  know  about  f to  improve 
our  sampling  procedure,  for  (20)  tells  us  that  the  error  will  be  decreased 
if  we  distort  p into  a new  distribution  with  more  mass  centered  in  the 
regions  where  f(x)  is  large.  In  practice,  it  is  not  hard  to  construct  a 
population  with  this  general  characteristic. 

A simple  illustration  may  clarify  this 
discussion.  Suppose  f(x)  is  a function  of  the 

JL 

type  pictured,  and  we  wish  to  find  f f(x)dx 

X) 

by  a sampling  procedure.  This  is  the  same  as 

"l  0 * X * 1 

finding  Eo(f(x))  where  f(x)  = ' $ 

0 otherwise 
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i.e.  we  pick  values  of  X from  a rectangular  distribution  on  (0,1),  and 
find  the  average  value  of  f(x).  It  is  clear  that  most  of  the  data  are 
wasted  in  this  procedure,  since  the  major  contribution  to  the  expectation 
comes  from  the  few  values  of  X which  we  happen  to  choose  from  the  neighbor- 
hood of  0.5,  One  easily  sees  that  a better  procedure  would  be  to  bias  the 
sampling  scheme  so  that  most  of  the  points  sampled  lie  in  the  vicinity  of 

0. 5.  It  is  then  necessary  to  weight  the  results  somehow,  so  as  to  compen- 
sate for  the  bias  introduced.  The  previous  paragraph  shows,  quantitatively, 
how  to  do  this.  The  term  "importance  sampling"  arises  from  the  fact  that 
the  amount  of  sampling  time  we  spend  in  studying  various  parts  of  the  popu- 
lation is  proportional  to  the  importance  of  those  parts  in  determining  the 
average . 

The  applicability  of  this  technique  to  our  problem  is  immediate.  In 

f | dS/r(t) 

the  case  of  the  hydrogen  atom,  we  are  interested  in  E \ e 

1. e.  we  sample  from  a space  whose  elements  are  paths  rft) , and  compute  the 

«f  ar/rCc) 

average  value  of  the  functional  F(r)  = e .In  practice,  the 

space  of  paths  from  which  we  sample  contains  only  a finite  number  of  ele- 
ments r^  (i  ® 1,  •••  , N),  since  the  paths  are  determined  by  the  choice  of 

a finite  number  (300)  of  Gaussian  deviates,  which  themselves  have  a basic 

—3  th 

increment  of  10  . If  p(r^)  is  the  probability  associated  with  the  i 


path,  we  can  write 
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N 

ill  F(r±)  P(ri)  = Ep  (F(r))  . 


We  have  observed  in  Section  IY  that  the  major  contribution  to  this  expecta- 
tion comes  from  a few  paths  r^  which  stay  near  the  origin,  i.e.  for  which 
F(r.)  is  very  large.  It  follows  from  the  foregoing  discussion  that  we  can 
decrease  the  error  in  our  estimate  by  associating  a new  set  of  probabilities 
q(r^)  with  the  paths.  The  q(r^)  will  be  an  improvement  over  the  p(r^)  if 
q(r^)  * p(r^)f(r^)  where  f(r^)  is,  roughly,  proportional  to  F(r^)  (it  is 
sufficient  that  f be  large  when  F is  large,  and  small  when  F is  small). 

Then,  as  before,  we  have 


(21)  Ep(F)  = -l  F(ri)p(r.)  = & [F(r.)  ^ ] q(r.)  = Sq(F^-)  , 

but  the  variance  of  the  q-scheme  is  less  than  that  of  the  p-scheme. 

There  arises  the  practical  question  of  how  to  find  a distribution  q 
having  the  desired  property,  without  hopelessly  complicating  the  calcula- 
tions. The  scheme  presented  here  is,  perhaps,  the  simplest  possible  one. 
Essentially,  we  bias  the  random  numbers  chosen  so  that  at  any  step  the 
’‘particle”  is  more  likely  to  step  tox^rard  the  origin  than  away  from  it. 
Most  of  the  paths  sampled  will,  therefore,  be  the  interesting  ones  x-jhich 
stay  near  the  origin.  Furthermore,  the  biasing  is  conducted  In  such  a 
fashion  that  the  factor  p(i\)/q(r^)  is  easily  computed.  ¥e  recall  from 
Section  IY  that  in  the  unbiased  x^alks  one  step  (say  the  k ) consisted  in 
choosing  three  numbers  from  the  Gaussian  deviates  with  mean  0, 


N 


P(r±) 
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variance  1,  and  basic  increment  .001.  These  three  numbers  represent  the 
components  of  the  step  on  each  of  the  Cartesian  axes . Then  the  probability 
p(r^)  to  be  associated  with  a 100-step  path  r^  determined  by  the  300 

choices  X1±,  Y1±9  Z1±,  •••  , xiQ0i*  Y100i*  Z100i  ls 


As  fan  as  we  are  concerned,  the  only  important  property  of  p(r. ) is  that  it 
has.  the  form 

100 

= Vi  * 10-9 


where  g is  a probability  density,  and  an  even  function.  Suppose  now  that 

we  bias  the  walk  in  the  following  fashion:  let  oC>/3>  o with  = 1; 

"fell 

at  the  k step,  pick  as  before  from  a Gaussian  distribution,  but  ignore 

the  sign  of  X,  ; then  we  play  some  game  of  chance  with  two  possible  outcomes 
K " k-1 

A and  B,  with  nrobabilities  <=^  and  /3 ; if  .ZL  X.  ^ o and  the  outcome  is  A, 

k-X  0 3 

interpret  X,  as  negative;  if  .21  X.  ^ o and  the  outcome  is  B,  interpret  X, 

K J = 0 - - 


r a i 

as  positive;  if  ,.2^  X..  < o and  the  outcome  is  j ^ k interpret 


Xk  as 


I positive  \ 

| negative  j # 


Sirailarly,  choose  Gaussian  deviates  and  play  the  auxiliary 


game  for  the  Y-  and  Z-  coordinates.  Clearly,  the  effect  of  this  scheme  is 
that,  at  a given  step,  the  absolute  value  of  a particular  coordinate  de- 
creases with  probability  oL  and  increases  with  probability  . Hence,  the 
walks  are  systematically  biased  toward  the  origin.  The  two-outcome  game 
involves  little  extra  effort;  if  we  take  = .6  and  /3  = .1;,  we  need 
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merely  pick  a digit  out  of  an  equidistributed  population  of  zeros , 
ones,  •••  , nines.  If  the  digit  chosen  is  not  greater  than  5,  we  call  the 
outcome  A,  otherwise  B. 

Let  us  write  down  q(r^) , the  probability  of  following  the  path  r^  under 
this  biased  system.  Let 


8 = 


k-i 

ft  if  has  the  same  sign  as  ^2^ 

k-1 

if  has  opposite  sign  to  j5l  x3 


(Similarly  9 (Yk)  , 9 (Zk))  . 

t 

It  is  easily  seen  then  that 


q(ri}  = [2^Xki>  9 t\i)]  9 (yid}]  [2«^>  9 (Zki}] 


Then 


/ \ 100  n OAA  -N.  „ N.  -300 

p(r. ) OAn  t | ■ 1 o-300  > •-  i /x  i 

iTJT7  = 2 31  ^ki59^9^  = ^ 


where  IL  is  the  number  of  times  the  auxiliary  game  yielded  the  outcome  A. 


We  can  rewrite  this  last  result  as 

p(r.)  \W_. 


pvriJ  /I  Vi  /_i\ 

^7  W v2  ft) 


X 300-I'k 


The  computation  of  this  factor  adds  little  complexity  to  the  IBM  work; 
we  merely  keep  a tally  on  N^,  the  number  of  outcomes  A,  and  read  the 
weight  factor  p(r^)/q(r^)  from  a table.  The  extra  IBM  time  involved  in 
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using this  importance  sampling  scheme  is  only  about  ±Q°/o  or  15%  . 

As  a preliminary  experiment,  mainly  to  test  the  qualitative  features 
of  the  method,  65  hydrogen  random  walks  with  importance  sampling  were  per- 
formed. ^ was  computed  by  means  of  equations  (10)  and  (21),  with 
t^  = 9 and  t^  = 16.  The  results  seem  quite  promising. 

One  point  concerning  our  treatment  of  the  data,  however,  should  be  men- 
tioned at  the  outset.  The  possibility  exists  that  we  may  sample  a path 
whose  weight  factor  p(r^)/q(r^)  is  exceptionally  large  because  the  auxiliary 
game  yielded  the  outcome  B improbably  often.  If  we  sample  enough  paths,  the 
effect  of  such  aberrations  will  be  negligible  (unless,  of  course,  the  im- 
portance scheme  was  incorrectly  chosen  and  virtually  neglects  significant 
regions  of  the  sample  space,  which  make  themselves  felt  only  through  such 
aberrations)  . However,  if  we  sample  only  a few  paths  and  happen  to  include 


one  with  a veiy  large  weight,  the  average  value  E 


{ft 


is  determined 


largely  by  the  single  path  with  high  weight,  and  the  rest  of  the  samples 
are  useless.  It  seems  reasonable,  in  a small  sample,  to  ignore  such  a path 
(certainly  so  if  we  know,  as  in  the  present  case,  that  it  represents  an 
unimportant  region  of  the  sample  space) . Consequently,  we  have  neglected 
a path  which  occurred  with  a weight  factor  25  times  greater  than  that  of 
ary  other  path,  and  100  times  greater  than  that  of  all  but  three  other  paths. 
These  other  three  paths  could  have  been  neglected,  with  almost  no  change  in 
our  results.  ("What  we  have  done  is  statistically  quite  justifiable,  and 
corresponds  to  the  rejection  of  measurements  which  fall  several  standard 
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deviations  from  the  mean!  the  measurements  may  be  correct,  but  the  sample 
is  not  large  enough  to  give  them  meaning.)  The  following  discussion  refers 
to  6b  random  walks , omitting  the  one  with  the  large  weight . 

On  the  basis  of  6U  hydrogen  random  walks  with  importance  sampling, 

« *5>l8  (true  value  ,5>00)  . We  recall  from  Section  IV  that  three  groups 
of  900  random  walks  without  importance  sampling  yielded  the  values 
^ = «30,  .36,  and  .29  respectively.  The  question  immediately  arises  of 

whether  the  accurate  result  yielded  by  importance  sampling  may  not  be  mere 
luck,  like  the  result  for  the  first  900  unbiased  walks.  To  Investigate  this 
possibility,  the  6I4.  walks  were  divided  into  arbitrary  subgroups,  and 
was  recomputed  on  the  basis  only  of  the  data  in  each  subgroup.  If  our 
result  were  due  to  luck,  the  computed  values  of  A-  should  vary  greatly  from 
one  subgroup  to  another.  Actually,  a high  degree  of  consistency  is  ob- 
served in  the  values  of  A^  so  obtained.  For  example,  a division  of  the 
walks  into  four  arbitrary  groups  of  16  yielded  the  values 


\ = 0.1*9  , 0.39  , 0.58  , 0.53 


respectively.  The  corresponding  computation  over  six  arbitrary  groups 
each  consisting  of  300  unbiased  random  walks,  yields 

X = 0.29  , 0.140  , 0.38  , 0.1*7  , 0.1*9  , 0.35 


respectively.  We  can  therefore  safely  state  that  groups  of  16  random- 

walks,  biased  according  to  our  scheme,  give  at  least  as  good  an  estimate 

of  (by  the  mean  value  method)  as  do  groups  of  300  unbiased  walks. 
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Importance  sampling  is  also  applicable  to  the  M slope  method11  of  Sec- 
tion V.  In  constructing  a histogram  for  unbiased  walks,  we  plotted  x 

versus  the  number  of  observed  paths  r.  for  which  x < F(r. ) - x + A.  If 

i x 

enough  paths  are  observed,  this  latter  quantity  is  proportional  to 

p(rj 

2 p(r.)  which  is  the  same  as  2>,.  . . —7- — r a(r.).  But  the 

x<F(r^)-x+A  x 1' 

latter  quantity  can  be  estimated  by  the  biased  walks,  and  is  in  fact  propor- 

4 

tional  to  the  sum  of  the  weights  p(r^)/q(ru)  of  the  paths  yielding  values 
of  F(r^)  in  the  range  (x,x  + A)*  Thus,  we  can  plot  a histogram  for  the  un- 
biased walks  on  the  basis  of  data  obtained  from  biased  walks ; for  a fixed 
number  of  walks,  the  biased  scheme  does  not  give  a good  picture  of  the  en- 
tire histogram,  but  yields  an  accurate  picture  of  the  tail,  where  most  of 
the  data  are  concentrated.  Such  a histogram  was  plotted  on  the  basis  of 
the  6H  biased  walks;  there  was  insufficient  data,  however,  to  define  any- 
thing more  than  the  order  of  magnitude  of  the  slope  of  a logarithmic  histo- 
gram, which  yielded 

With  very  little  extra  computation,  the  biased  hydrogen  walks  furnish 
us  with  data  for  biased  haimonic  oscillator  random  walks . We  recall  that 
the  lowest  eigenvalue  of  the  harmonic  oscillator  is  computed  from  the  cuan- 

-J,fcx2(t)dtr 


tity  E > e ° 


* , which  we  estimate  by  random  walks  • The  walks 


which  spend  much  time  near  the  origin  contribute  the  larger  values  of 
this  exponential;  hence,  the  biasing  scheme  for  hydrogen,  which  stresses 
these  same  paths,  might  be  expected  to  be  useful  for  the  harmonic  oscilla- 
tor. On  the  basis  of  100  biased  harmonic  oscillator  random  walks,  we 


■ 


■ 


■ 


■ 

. '{■•..  •"  ■ • : - ■ • 5 

■ 

. ’ ••  ■ 


,*  ••  • > 

. . 


. 

V ••• 


. 


•:  ■ 

. 


....  . 


■ ■ • . 
’ ’ • 


. 


...  r l . 


obtain  * 0.81  (true  value  0.71).  Division  of  the  data  into  ten  sub- 


groups,  each  of  ten  walks, 

yields 

■ 0.66 

0.85 

0.86 

0.78 

0.66 

0.8? 

1.15 

0.96 

0.81 

0*71 

respectively,  for  each  of  the  subgroups*  In  this  case,  our  final  result 
for  is  less  accurate  than  that  obtained  from  unbiased  walks,  but  our 
subgroups  are  more  consistent. 

One  can  think  of  many  other  types  of  importance  sampling,  e.g,  where 
the  strength  of  the  biasing  is  a function  of  position,  or  where  the  biasing 
scheme  is  successively  improved  by  some  sequential  procedure.  Even  with 
the  present  biasing  scheme,  lowest  eigenvalue  of  helium  should  be  obtainable 
to  within  2rfo  in  two  hours1  time  on  the  SWAC.  The  authors  will  welcome 
correspondence  concerning  unsolved  potentials,  to  which  the  sampling  method 
may  be  applicable* 


The  authors  wish  to  express  their  thanks  to  Dr.  Wolfgang  Wasow,  of 
the  National  Bureau  of  Standards,  whose  valuable  comments  brought  to  light 
many  obscurities  in  the  original  manuscript;  and  to  Dr.  Everett  C.  Yowell, 
of  the  National  Bureau  of  Standards,  who  not  only  supervised  the  IBM  work 
infallibly,  but  reduced  card  jams  as  fast  as  the  junior  author  could 
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APFE1DIX  A 

COMPUTATION  0?  DISCRETIZATION  ERROR  FOR  THE  HARMONIC  OSCILLATOR 


The  discretization  error  K 0 was  defined  in  Section  III  as 

r 

nt 


n 


< exp  j - \ j£  S,J  - <^exp  | - f x2(z)dV 


We  define 


<Zej®  | - J*  x2(E)dT  | 

X(n,t)  - <<xpj-  \ Sk2}>  . 

We  shall  compute  l(n,t)  exactly.  Since,  by  equation  (5a) 

lira  I(n,t)  = /exp  <f  ~ J x2(r)dr  V/>  , 

n— >oo  \ l o J ' 

our  calculation  of  l(n,t)  will  yield  the  exact  value  of  the  Wiener  integral 

as  a byproduct.  Then  we  can  compute  K 0, 

x 

We  recall  that  S,  - S,  n is  a Gaussian  random  variable  with  mean  0 and 

k k-1 

variance  1 . Therefore 

Prob{xl  < Sx  < ^ + d^,  x,  < S2  < X2  + a**,  •••  , xnt<  Snt  < Xnt  + dxntj 


nt 

2 


(2i»)  exp 


{-![ 


*1  + + (^-X2}  + *"  + 


0} 


dx_  dx^,  • a ’ dx  , 

x 2 nt 


Consequently, 


' 

- ■ ■ > 


• . ..  ■; 


■ 


. 
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I(n,t)  - /exp  { - \ X Sfc2  j 

f 1 2 ir  2 2 

exn  - — kJi  \ " 2 L^l  +(x2~xi)  + ••• + 
l n 


dxi  dxnt  . 


(x  , “X  , i ) 
nt  nt-1 


Integrals  of  the  form 


I = (TT) 


m 

2 


m 


Z a.  . x.x. 

, 1,0=1  10  1 0 


dx^ 


dx 


m 


(where  a.  . is  symmetric  and  positive  definite)  are  easily  evaluated  by 
10 

transformation  to  principal  axes.  By  a pure  rotation  of  the  form 

- det  B * 1)  we  can  bring  I into  the  form 


y « B x (where 


The  A’s  are  the  elements  of  the  diagonal  matrix  A - B A B~*^. 
Hence 


I = (det  A ) 


-1/2 


(det  (a.  .) 

V 10 


Applying  this  result  to  l(n,t),  we  obtain 


Let  us  consider  the  k x k determinant 


Expanding  according  to  the  first  row  or  column,  we  establish  the  recur- 
sion relation 

Dk  (a,b)  • a Dj{_1  (a,b)  - (a,b)  . 

If  we  define  D * 1 and  D-  * b,  this  equation  is  true  for  k * 2* 
o x 

We  rewrite  the  recursion  relation  as 


Vl  + 


k-1 


a D, 


(k  ^ 1) 


■ 

. ' ’ 


■ 


- 

■ 

: I 

« 

• 

. 

■* 
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Utilizing  the  identities 

cosh  (k  + 1)  9 + cosh  (k  - l)  0=2  cosh  9 cosh  k 9 
sinh  (k  + l)  0 + sinh  (k  - l)  9 = 2 cosh  9 sinh  k 9 
we  can  write 

Dj  (a,b)  = A cosh  k 9 + B sinh  k 9 


where 


cosh  9 = a/2 


To  fit  the  starting  conditions  on  D and  Dn  we  set 

o 1 


A = 1 B = 


Thus 


I( 


n,«  - [ 


cosh'  (n  t 9)  + 


2b- a 
/ a2- It 


1 


sinh  (n  t 9)  J 


-1/2 


where  cosh  9=1+  1/n  * Solving  for  e , ve  obtain 


0 


= 1 + l/n2  + i \fz  + l/'n2  = l + 'i--+\+0  (n-3) 

1 n 


As  n — =>  oo  , e 


n t 0 


e"^^  « Therefore 


lim  I(n5t)  = <^exp  - J x^CEjd'cJ’^  = £cosh  t yfe 


' 


* 


• • , • * ■ 
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In  estimating  K 2>  we  shall  coniine  ourselves  to  the  case  t > > 1, 
x 

t/n  <<  1.  This  is  the  case  of  interest  for  the  harmonic  oscillator,  where 


we  took  t = 5,  n = 100.  Therefore  we  neglect  terras  involving  the  factor 

-t\/2 


e ^ compared  with  terras  involving  e"^v  ^ . Thus 


• /yp  j “ I x2(T)dtjy> 


v'T 


l(n,t)^  e 


ntQ 

2 


rx* — i — i 
L 


-1/2 


-W2 


We  shall  expand  the  difference  e ^ in  powers  of  l/n,  retaining 

the  term  of  lowest  order.  It  is  convenient  to  deal  with  logarithms  at 


this  point 


The  vanishing  of  terms  of  order  n in  the  expression  for  9 is  signifi- 
cant, since  it  implies  that  the  expression  for  log  l(n,t)  contains  no 
error  term  of  order  t/n.  This  is  the  computational  reflection  of  our 
qualitative  conjecture  that  the  error,  to  first  order,  involves  only  a 
shifting  of  the  time  origin,  rather  than  a dilation  of  the  time  scale. 


. 


. 


' 


■ 


■ 
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We  obtain,  then 


log 

I(n,t)  _ -VS 

_ JL.  + 

^ 2 

2\/Sn 

-WS 

I(n, 

.*)  - s”^  Ti 

— - ^ 4. 

L 

2\/5n 

Hence 


K p = -i—  + 0(t/n' 
x^  2\/Sn 


0(t/n2)  , 

0(t/n2)l  \/2 

) . 
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APPENDIX  B 

FURTHER  STATISTICAL  PROPERTIES  OF  THE  RANDOM  WALKS 


Since  many  properties  of  the  hydrogen  random  walks  are  predictable 
analytically,  it  seemed  worthwhile  to  check  the  experimental  results  against 
some  of  these  predictions. 

The  most  significant  result  obtained,  and  one  which  may  represent  a 
source  of  error,  concerned  the  distribution  of  the  end  points  r(l)  of  the 
hydrogen  walks.  We  recall 


//Too  T2  TToo  TS  7Too  \7 
■<«  “ + {v* iz’-> 


( 


100 


Since  \k»l  XjJ  , etc,,  are  Gaussian  with  mean  0 and  variance  100,  then 

2 2 2 2 
[r(l)r  has  the  same  distribution  as  X + Y + Z where  X,  Y,  and  Z are  in- 
dependent Gaussian  random  variables  with  mean  0 and  variance  1. 


Hence 


2 2 2 
z 


n p r ~ x +y  * 

Prob  / r(l)  <^-  1 = (2lT)  J J J e 

l J 2 2 2 2 


2 2 2 2 
X +y  +z 


dx  dy  dz 


= (27 JT3^2  f e 2 U1T  r2  dr 
o 


2 oc  - c , » 

= erf  - — • — e = F(oc) 

^tr 


where  the  last  line  is  obtained  through  an  integration  by  parts. 


; 

. 


• r "•••  • . ~ : - ■’ 

, 

••  ' ^ ■ . 


* 


■ 


- 
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If  we  tabulate  the  end  points  r(l)  of  the  first  900  random  walks,  the 
expected  number  of  walks  with  r(l)  4s*  is  900  F(=*l) . Let  us  define 

S(°0  * number  of  walks  with  r(l)  (empirical) 


It  Is  easy  to  show  that 


<[■ 


]> 


s(oc)  - 900  F(=*-)  I > s (7-  (ex.)  = 900  F(o^)(l  - F(cX-)) 

r ®<-)  "I 

Furthermore,  the  probability  that  |s(®*)  - 900F(C<-)  | 2c7_('>c)  J-  is 

3=K<^)J 


r*68 1 

is  4 .95  I . 

( .997 J 


Hence,  discrepancies  of  more  than  2 cr(c*> ) between  the  empirical  and  theor- 
etical distributions  are  to  be  regarded  with  suspicion.  We  compare  the  two 
distributions  in  the  table  below* 


DISTRIBUTION  OF  END  POINTS  OF  900  HYDROGEN  WALKS 


F(c=c) 

900  F(<*) 

s(<*0  • 0 

• (pc) 

,S(^)-900  F(«-)  , 

1 {<=< ) 

o*5 

.030 

27 

28 

5.1 

0.2 

1*0 

.197 

177 

185 

11.9 

0.7 

1.5 

.U78 

130 

U27 

i3.o 

0.2 

2.0 

.739 

665 

6U5 

13 .2 

1.5 

2.5 

.900 

810 

753 

9.0 

6.3 

3.0 

.971 

87U 

836 

5.0 

7.6  • 

3.5 

,99k 

89U 

867 

2.3 

11.7 

U.o 

.999 

899 

880 

1.0 

19 .0 

The  tail  of 

the  empirical 

distribution  S( 

evidently  falls  off 

much  less  rapidly 

than  that  of 

a true  Gaussian, 

, The  magnitude  of  the 

. • *'  . ■ r . • . , " ' -•  . • V 

-v  . : 'i  fi  i:  f,  'V  i ' o »<flc 


■ 


i 


■ • 
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• 
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\ 

.... 
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■ 
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resultant  error  in  and  possible  remedies,  have  not  been  investigated, 

mainly  because  the  SVJAC  will  not  employ  the  Rand  Corporation  Gaussian 
deviates.  Qualitatively,  the  effect  of  the  no  n-Gaus  si  an  behavior  is  to  in- 
crease the  values  of  r(£),  and  thus  to  decrease  the  values  of  J dt/rfC)  * 

- • o 

This  may  provide  a partial  explanation  of  another  discrepancy  between  theory 

and  observation,  discussed  below. 

Comparisons  were  also  made  between  the  theoretical  and  experimental 
mean  values  of  integrals  arising  in  the  hydrogen  and  helium  computations. 
Since  the  operations  of  averaging  and  summing  (or  integrating)  commute,  we 
have 


|x  at/r(-c)\>  - 


But 


r(x) 


- (2fTX)“3/2  f e-r2/2r  U7T  r2 


, /2  /l 

r v IT  V r 


Hence 


<f  **«>  - ( * 


-1/2 


dt 


x/Trr 


= 1.59 


The  empirical  value  for  1800  samples  was  1.U3. 


In  order  to  appreciate  the  significance,  if  any,  of  this  discrepancy, 

2 p1 

one  must  know  the  variance  O**  of  a single  observation  of  I dT/rfC) . 

r,  O 

Knowing  , we  then  know  the  variance  of  the  average  of  1800  observa- 

f ■ ° ^2  ^o2 

tions , since  o*  = . t^qq • 


- 

' 

. r 

■ 

. 


. 





4 


3 y definition 


^0 


/,! 
W 

Wo 


3-' 2 - ( P d't/r(r)  - — i- 


v2  nr 


)2>-  <(j>/W2>-^ 


a 


If  we  think  of  f dX/r(?Z)  as  a Riemann  sum  of  many  non-independent  random 


± 


variables  , it  is  easily  seen  that 


r?Eh; 

M 


\ 


d :/r\Z) 


/ (2tr)3  J0  Jc  ?372  d\  dr2  * 


CD 


f 6 r. ^ 

U rjz 


*1  +7l+*l 

2T1 


2 2 

(xp-x,)  +(y9-y1)  +(s2-z1) 

2<Vv 


2 2 >n  2 2 


•00  /x,  +y  +Z-,  /x0  +yn  + z9 

V a.  I.  I V ^ £ 


dx,  •••  ds0 

J.  u 

= U log  21 


Thus 


<yf  = I*  log  2 - JL  = .22 
cr0 

<r  = - — 2~~  - .on  . 

\TTHoo 


(In  good  agreement  with  the  empirically 

o 

observed  <T^ " - .19*) 

0 


^Evaluation  of  th$ 


integral  is  left  as  an  exercise  for 


the  reader. 


« 


- 56  - 


The  difference  of  *16  between  the  calculated  ana  ■ .. crved  value  of 

■i 

J d'C/rff  ) is  therefore  equal  to  l5  which  is  far  t eyond  the  limits  of 

statistical  inaccuracy.  Dr.  R*  P.  Feynman  first  called  attention  to  a 

major  systematic  error,  and  gave  the  following  method  for  estimating  ite 

r1  , , . 

What  we  actually  sample  is  not  the  integral  j d /r ('■'),  but  rather  the 


10( 


u 

o 


Riemann  sum  *01  ^ — — 

^loo5 


Thus,  the  average  to  be  considered  is  that  of 


the  Riemann  sum,  which  is  somewhat  smaller  than  that  of  the  integral  (the 
error  encountered  here  is  the  ’'discretization  error*1  discussed  in  Section 


It  is  geometrically  evident  that  a fairly  good  estimate  of  the  Riemann  sum 
in  parentheses  is 


i 


1.005’ 

oo5 


ar 

v'T 


1.003  - .071 


1.86 


Thus  <^.01  kJj_  ?'(Y7iooy  ^ = (compared  with  1.5#  for  <^  | df/r(f^>) 

and  our  observed  value  1*U3  is  in  error  by  5 cr,  instead  of  1 The  re- 
maining error  may,  perhaps,  be  accounted  for  by  the  non-Gaussian  behavior 

of  the  tail  of  the  distribution  of  rfl)  (discussed  above),  and  the  conse- 

1 

quent  depression  of  values  of  f d'^/r ) . 


’ 
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One  can  also  calculate  the  average  value  of  the  integral  of  the  helium 
integration  term. 


1 


•\f  (xg-xp  2+  (^-yp  ^(Zn-Zj)2 


2 2 2 

*1  *?1  *Z1 
2X 


>y,  +z, 


2t 


dx-j  • • • 


da, 


The  observed  value  was  .9 3.  The  discrepancy  is  probably  to  be  explained 
along  the  same  lines  as  above. 


APPENDIX  C 


STATISTICAL  METHOD  FOR  FINDING  THE  LOWEST  EIGENFUNCTION 


The  reasoning  in  Section  II  leads  easily  to  a statistical  method  for 
finding  (x) , the  lowest  eigenfunction.  From  Section  II  we  have 

~ A.  ,t 

P(x, t)dx  - H C-  .(x)e  J dx  . 

0 J 

Furthermore 


P(x,t)  cte  « 

Prob  i particle  survives  till  time  t and  is  in  (x,x  + dx)  at  time  t r 

1 r i J 

Prob  'I  particle  picks  a path  such  that  x < x(t)  < x + dx  j x 
Prob  ^ particle  survives  till  t,  given  that  x < x(t)  < x + dx  J 


But 


Prob  I particle  picks  a part  such  that 


. . It 

x < x(t)  < x + dx  l 


v2Tf  t 


-x  /2t  . 
e ox  $ 


and  the  probability  of  survival  till  t,  given  that  x < x(t)  < x + dx,  is 
just  the  average  of  the  survival  probability  exp  ^ / V(x(£)dTj^  over 


the  set  of  paths  such  that  x < x(t)  < x + dx.  Hence 


P(x,t)dx  = 


1 -x2/2t  x 

— ..  A ' / 


3x  <(eoq>\  -J  V(x(t))dZ 


x < x(  t)  < x + dx 


> 


V2  Tf  t 

The  conditional,  average  in  brackets  can  be  found  empirically  by  subdividing 
the  x-axis  into  short  intervals,  and  computing  the  average  of 


' 


. 

■ 

• • 


. 


' 


P(x,t)dx  for  large  t,  we  can  thus  determine  b$  (up  to  a normalizing 
factor)  . By  considering  several  different  values  of  t and  waiting  until 
the  rate  of  time  decay  of  P[x9i)  becomes  the  same  at  all  points  x,  we  can 
ascertain  when  the  space  distribution  has  settled  down  to  the  correct  eigen- 
function. This  settling-down  process  can  be  speeded  by  starting  the  walks 
not  from  a delta-function  distribution,  but  rather  from  a distribution  which 
approximates  our  best  guess  for  (x)  . Then  the  coefficients  0o,  ••• 

will  be  small  • 

No  experiments  have  yet  been  performed  with  the  above  scheme,  which 
clearly  requires  very  many  random  walks . 

If  we  start  our  walks  from  a point  x^  instead  of  the  origin,  then 

,t 

P(x,t)  = 2 ) Y.(x)e  3 

o u J 


and  hence 


N^(xo)e  = j P(x,t)  ^f/1(x)dx  = <Qaxp  j - V(x(x))dT  J*  v|/^(x(t)^ 

If  we  have  an  empirical  estimate  of  and  An  from  an  earlier  set  of 

walks,  we  may  perform  a new  set  of  walks  starting  at  x and  use  the  last 


xIf  N paths  x^(tr)(i  =1,  •••  , N)  are  sampled,  an  equivalent  and  com- 
putationally easier  estimate  of  P is 

P(x,t)dx  = pT  2 exp*  { - I V(x.$))dt \ . 

Nx<x,(tVx+dx  L % x J 


. 

, 


. 

jjpir 
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for  (x, ) . The  osef nines s of  th: 


pro' 


equation  to  improve  our  gueg^  a.vx  5 - . 

cedars  is  doubtful,  since  it  depends  on  the  ,!  approximate  orthogonality”  of 
our  final  ‘^le  true  ''V2>  "¥y  etc.,  and  involves  an  inordinate 

amount  of  effort  to  improve  our  guess  for  at  a single  point* 
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